Unconventional superconductors under rotating magnetic field II: 

thermal transport 
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We present a microscopic approach to the calculations of thermal conductivity in unconventional 
superconductors for a wide range of temperatures and magnetic fields. Our work employs the 
non-equilibrium Keldysh formulation of the quasiclassical theory. We solve the transport equations 
using a variation of the Brandt-Pesch-Tewordt (BPT) method, that accounts for the quasiparticle 
scattering on vortices. We focus on the dependence of the thermal conductivity on the direction 
of the field with the respect to the nodes of the order parameter, and discuss it in the context 
of experiments aiming to determine the shape of the gap from such anisotropy measurements. 
We consider quasi-two dimensional Fermi surfaces with vertical line nodes and use our analysis to 
establish the location of gap nodes in heavy fermion CeCoIns and organic superconductor k-(BEDT- 
TTF)2Cu(NCS)2. 
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I. INTRODUCTION 

In the preceding paper—, hereafter referred to as I, we 
developed a theoretical approach to the vortex state in 
unconventional superconductors that allowed us to ob- 
tain a closed form solution for the equilibrium Green's 
function, and therefore efficiently compute the density 
of states and the specific heat for an arbitrary orienta- 
tion of the magnetic field. In this work we extend our 
approach to the calculation of transport properties, de- 
velop the formalism for computing the electronic thermal 
conductivity, and compare our results with experiment. 

The rationale for both calculations is to provide the- 
oretical guidance and support to continued attempts to 
establish the measurements of the anisotropy of the spe- 
cific heat and thermal conductivity under rotated field as 
a leading tool in determining the structure of the energy 
gap in unconventional superconductors. While a num- 
ber of techniques test the symmetry of the gap via the 
surface measurements, the corresponding bulk probes are 
few. Semiclassical treatment of the quasiparticle energy 
in the vortex state incorporated the Doppler shift due to 
local value of superfiuid velocity associated with the cir- 
culating supercurrents. This approach predicted that, at 
low fields, the density of field- induced states at the Fermi 
surface oscillates as a function of the field direction and 
has a minimum when the applied field is aligned with the 
nodal direction, |A(p)| = 0, hence the suggestion to use 
the measurements of the low temperature specific heat to 
determine the position of nodes^i. The experiments are 
quite challenging, and, for now, have been carried out in 
few matcrials^i^i^. 

Variations in the density of states also influence trans- 
port properties, and the measurements of the electronic 
thermal conductivity under a rotated fleld have been used 
more extensively to study unconventional superconduc- 
tors and infer the gap structure^'^'^ d^d^d^i^'^d^ . Exper- 
imentally, for a fixed direction of the heat current and 



rotated field, the dominant twofold anisotropy is that 
between the transport normal to and parallel to the vor- 
tices; a much smaller signal is attributed to the existence 
of the nodes, see Ref. [Tj for a recent review. Theoretical 
analysis of the thermal conductivity is also much more 
challenging. There are conceptual difficulties with ex- 
tending the "local" semiclassical approach to calculations 
of the response functions, especially for clean systems 
where the mean free path exceeds the typical length scale 
for the variations of the superfiuid velocity, the intervor- 
tex distance. Even in moderately dirty systems, where 
the use of semiclassical method is justified, it yields, at 
best, a local value of the thermal conductivity, which 
varies from point to point; consequently the averaging 
procedure to obtain the experimentally measured value is 
far from obvious. Naive averaging completely misses the 
twofold anisotropjJ^, and therefore is not trustworthy. 
The semiclassical approach does not naturally include the 
scattering on the vortices, and attempts to introduce it 
phenomenologicall y^^i are promising, but have not yet 
lead to a consistent description. Moreover, the exper- 
iments on all but high-Tc and some organic supercon- 
ductors are done at fields that are a significant fraction 
of the upper critical field, Hc2, where the accuracy of 
the semiclassical approximation may be called into ques- 
tion. Consequently, we argued that a more microscopic 
approach is needed. 

We use a quasiclassical version of the Brandt-Pesch- 
Tcwordt approximatio n^^! , where the normal electron 
part of the matrix Green's function is replaced by its spa- 
tial average over a unit cell of the vortex lattice. Remark- 
ably, this approximation allows for the closed form solu- 
tion for the Green's function that we found in Refs.[lll20l. 
and used, with a fully self-consistent treatment of the or- 
der parameter and impurity scattering, to determine the 
behavior of the specific heat across the T-H phase dia- 
gram. Below we review these results and develop a linear 
response theory for thermal transport. Implementation 
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FIG. 1: (Color online) The model considered in this paper. 
Calculation of the thermal conductivity is done for a quasi- 
cylindrical Fermi surface, when the heat current or tempera- 
ture gradient and (rotated) magnetic field are in the afc-plane. 
The order parameter is assumed to have a d-wave symmetry. 



of the approximation in the framework of transport-like 
quasiclasical (Eilenberger) equations^ii^ ensures that we 
account for the difference between single particle and 
transport lifetimes in scattering off of vortices: the char- 
acteristic intcrvortcx distance is large compared to lattice 
spacing, and the scattering on the vortices corresponds 
to small momentum transfer, and hence forward scat- 
tering is important. This allows us to treat the twofold 
anisotropy (transport normal and parallel to the vortices) 
on equal footing with the effect of the nodes, and to de- 
velop a consistent picture of the behavior of the thermal 
conductivity and the specific heat under the same as- 
sumptions. 

As in part f, wc consider a quasi-two dimensional 
Fermi surface to focus on the comparison with the data 
on heavy fermion CeCoIns. fn that materials the spe- 
cific heat data were interpreted as supporting the d^y 
gap symmetry^; however, as pointed out in our short 
communication^^ and preceding paper, the anisotropic 
part of the specific heat changes sign at moderate fields 
and temperatures, with maxima, rather than minima, 
for the field aligned with the nodes. At low T and H, 
in the region of validity of the semiclassical method, our 
results agree with that of the calculations utilizing the 
Doppler shift approach, with minima for the field along 
the nodes. Consequently, in light of these observations, 
we reinterpreted the results of Ref . as possibly support- 
ing the dx'i-y'^ gap symmetry. Similar gap structure was 
inferred by fzawa et al. using phcnomcnological interpre- 
tation of the thermal conductivity measurements^. We 
provide the detailed analysis of the thermal conductivity 
here. 

In the following section we briefly review the approach 
and the main results of the preceding paper, such as the 
closed form expressions for the Green's functions nec- 
essary for computing linear response to the gradient of 
temperature. Sec. IIIII gives the derivation of the ther- 
mal conductivity using Kcldysh formulation for the non- 



equilibrium theory of superconductivity, with some de- 
tails relegated to Appendix VK\ As in I, we find that the 
simple example of a 2D d-wave superconductor with a 
cylindrical Fermi surface provides a semi-analytically ac- 
cessible path towards understanding some of the crucial 
features of our results, and consider it in Sec. IIVI Sec. |V] 
is devoted to the calculations for more realistic quasi- 
cylindrical Fermi surface (Fig. [T|), and at the end of it 
we discuss the results, compare them with the data on 
CeCoIng and organic k - [BEDT) - TTF. 

Our Sees. IIVI IVl and I VII are intended for those readers 
who are interested only in the overall physical picture 
and the behavior of the measured properties; the figures 
in Sec. |V] show the main differences between the self- 
consistent and non-self-consistent calculations. Finally, 
our conclusions provide a side-by-side comparison of the 
specific heat discussed in I and the thermal conductivity 
results, and outline implications for future experiments. 



II. QUASICLASSICAL APPROACH AND THE 
EQUILIBRIUM GREEN'S FUNCTION 

A. Basic equations and formulation 

We follow Ref. I in considering the quasiclassi- 
cal (integrated over the quasiparticle band energy) 
Green's function in a singlet superconductor in mag- 
netic fiel d^-^'^^1^^'^'^1^^'^^ . In the spin and particle-hole 
(Nambu) space the matrix propagator depends on the 
direction at the Fermi surface (FS), p, and the center of 
mass coordinate, R, and is 

5(H,fte,= (.4-^)^ (1) 

We write the quasiclassical equation for the real- 
energy, e, retarded, advanced, and Keldysh propagators, 
which enables us to carry out non-equilibrium calcula- 
tions, see Appendix S below and Refs. [Ull^Ili. The 
retarded (R) and advanced (A) functions g = g^'^ sat- 
isfy (we take convention e < 0) 

[(e-f ^v/(p)A(R))f3 - A(R,p) -5,™p(R;£), 

5(R, p; e)] + ivf (p) • g(R, p; e) = , (2) 

together with the normalization condition 

ff«'^(R,p;£„,)2 = -^2^. (3) 

Here v/(p) is the Fermi velocity at a point p on the FS. 
The vector potential A(R) describes the applied mag- 
netic field, and the self-energy, a, (different for the re- 
tarded and the advanced components) is due to impurity 
scattering. The mean field singlet order parameter, 

A-f ^'^^A^i 
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is self-consistently determined using the Keldysh func- 
tions 

A(R, p) = I ^ ( V{p, p') f''{R, p'; e)U • (5) 
In Eq.Q we used a shorthand notation 



(6) 



where n/(p) = Nf{p)/J\ff, with Nf{p) the density of 
states (DOS) at a point p on the Fermi surface in the 
normal state, and Mf = J o?Pfs^/(p) the net density of 
states. 

Throughout our work we use a separable pairing inter- 
actions, 



vip,p') = v,yip)yip'), 



(7) 



where 3^(p) is the normalized basis function for the angu- 
lar momentum representation, (3^(p)^)fs = 1- Hence the 
order parameter is A(R, p) = A(R)3^(p). For example, 
for a 2D dx2_y2 gap we take ^(p) — V^ipi —pl)- 

We include the impurity scattering via the self- 
consistent t-matrix approximation, with the self-energy 



(R;e) 



D + 11 ilJ2^imp 



(8) 

Here riimp is the impurity concentration, and, if u is the 
isotropic single impurity potential, the t-matrix is defined 



i{R.; e)^ui + uJVf {g(R, p; e)) fs ^ (R; e) . (9) 

The component labeled D 1 drops out of the equations 
for the retarded and the advanced Green's functions since 
the unit matrix commutes with the Green's function in 
Eq. This term, however, appears in the Keldysh 

part, and affects the transport propertie a^^'^° (see also 
Appendix |^ . Below we parameterize the scattering by 
the "bare" scattering rate, F = 
shift 6o of the impurity scattering, tan Sq 



a/irAff, and the phase 



B. Equilibrium Green's function 

In I we solved the quasiclassical equations in the vortex 
state. We took the superconducting order parameter in 
the form A{R, p) A{R)y{p), where 



A(i?) =^A„(R| 



(10a) 



(R|n)=^C(';)^-_$„(z,fc,). (10b) 

We showed that for J2ky l^fc"^l^ ~ ^ coefficients 
C^"'* determine the shape of the lattice, while A„ is the 



amplitude of the order parameter in the n-th Landau 
level channel. The expansion of A(i2) is in the Landau 
level function of the renormalized coordinates. 



$„(X, ky 



(11) 



The anisotropy factor, Sf, for the field applied at an angle 
9h to the z-axis, is 



Sf = J cos2 eH + -^ sin^ Oh 



(12) 



where v^^ = 2{y^{p)v'j^^{p^))^s and v^^^ = 
2(3^2(p)ujl (P2))fs; here W|| is the projection of the 
Fermi velocity on the z axis, and v±i with i = x,y is the 
projection on the axes in the plane normal to z. 

Following the BPT procedure we replaced the diagonal 
part of the Green's function, g, with its spatial average, 
and introduced the ladder operators for the Landau levels 
which allowed us to solve for the anomalous (Gorkov) 
hmctions in terms of g. 



/(R,p;£) = ^/„(p,e)(R| 



(13a) 



/m(p, e) = ig ^(-5_(p))™-" !?„,„(£, |p|)A„(p; e) . 

n 

(13b) 



Here 



■5±(p) 



Vf{p)^±ivf{p)y 



\vf\ 



(14) 



with 



Vf{p)x=Vf{p)J^/Sf , Vf{p)y ^Vf{p)y^ , 



(15) 



and 



1^^/ (P)l = ^Jvfip)l+ifip)l (16) 
The coefficients of the expansion arc given by 

min{m,n) / 2f A \ 

p„,„(e, ipi) = v^— y: (-1)'" , 



(17) 



withni(j) = j + (|m-n|-(m-n))/2, ?i2(j) = j + (|m- 
n\ + (to — n))/2 in each term and 



■ \ n 1+712 



V2 



(?! — ?ii)!ni!?i2! 



(18) 



where iy'^")(z) is the n-th derivative of the function 
W{z) = exp(— z2)erfc(— iz). 

We then use the normalization condition, 



(19) 
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in the spatially averaged form, with /i /2 the spatial aver- 
age of the product, to find the equilibrium Green's func- 
tion 



g = -iTr/vT+P, 
2 



(20a) 



fc.i>0 



X ( n |a'"a'| m ) I ^ 1 T/j/C^+'+i) 



k+l 



-1) / V2e 



, (20b) 



theory. We first derive the general formula for the heat 
conductivity tensor, based on the closed-form solution 
for the quasiclassical retarded and advanced propaga- 
tors found above, and using the non-equilibrium Keldysh 
approach. The Keldysh part of the full quasiclassical 
Green's function carries information about both the spec- 
trum and the distribution of quasiparticles, and the heat 
current is defined as energy transfer by quasiparticle o^^i^^ 



where w — ju^l/v^A, and the prime over the k,l- 
sum denotes the restriction that the matrix element 
{n ja^'^a'l m) = y/nlml/{n — k)l{m — 1)1 is non-zero only 
for k < n, I < m and k — I = n — m. This expression 
reduces to the form of g obtained previously if we trun- 
cate the order parameter expansion at the lowest Landau 

jj,,,j, ^19.25.31.32.33 



9 ■ 



(21) 



This latter form is useful for semi-analytical calculations. 



III. 



HEAT CONDUCTIVITY 



A. Linear response and thermal conductivity 

We now proceed to derivation and analysis of expres- 
sion for the thermal conductivity in the linear response 



j„(R) = 2iV/ J dp rif ip)vf{p) J ^£5^(p,R;£), 

— OO 

(22) 

where g^ is the diagonal component of the Keldysh prop- 
agator. 



In equilibrium j/j = as expected, and in linear 
response we define the heat conductivity tensor k via 
jh = —kVT. We linearize the equations to find the first 
order corrections to the retarded, advanced and Keldysh 
propagators, S'g^'"^'^ , with respect to VT. This implic- 
itly assumes that the inhomogeneity due to the tempera- 
ture gradient is much smaller and occurs on much longer 
scales than the inhomogeneity due to vortices, impuri- 
ties, etc., which is the case experimentally. The details 
of the derivation of k are presented in Appendix [X] and 
here we give only the final expression, 



T 

where we defined. 



Nf 


f de 




2 J 


' T T^cosh\e/2T) 


Gi = 


- 


1 


^ 2{g^ -g^) 


G2 = 






2(5«-5-4) 


G3 = 






2{gR-g^) 


G4 = 


_jja _ 


1 


^ 2{gR-g^) 



/%G 1 
dpFsn/(p) i'/,j(p)i'/j(p) 777; Fr7^( ^"^5^) (23a) 



[-(/« + /^)(a'' - A^) + (/«■ + /■^)(A«' - A^)] (23b) 



[(/«■ - + {t-t)iA^ - A^)] (23c) 

[(/«■ + !^){^ + A^) + [f + /^) (A« + A^)] (23d) 
[-(/«■ - /^)(a'' + A^) + (/«■ - /^)(A« + A^)] , (23e) 



and used the following notations: D°-(e) ~ D^{e) — We can re- write Eq. (|23ap as 

D'^ie), £"^(2) = S^(e) - I]'4(e), and A^'^(R,p;e) = 



A(R, p) -I- A^p(R;e). In both Born and unitarity scat- Kij{T,H) f de 2 £ 



+00 



■cosh-^— (24) 



tering limits D'^ = 0, which simplifies this resuHi^i^. T J 2T T'^ ''^'^^^ 2T 

— OO 

X J dpFsW/,,w/j N{T,H;p,e) th{T,H;p,£) 
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Here N(T,H;p,e) = iV/n/(p)(-Jm g^/ir) is the angle- 
dependent DOS, and th = —iG2/ (6*104 — G2G3) has the 
meaning of the transport hfetimc due to both impurity 
and vortex scattering. In the normal state th = Tn = 
1/27 (7 = Tsin^ So) and -3m g^/ir = 1. Notice that the 
transport and the single-particle lifetimes are different. 

Several limiting cases are useful for developing a quali- 
tative understanding of the physical picture. In the Born 
or unitary limit = = 0. If we truncate the ex- 
pansion of the vortex state at the lowest Landau level 
function, ?i = 0, and neglect the off-diagonal impurity 
self-energy Ai„ip = 0, we obtain from Eqs. (|13ll^ 
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= (25) 



. « 2V^A 
^9 -T^W 

\vf\ 




|Ao3^(p)|^ 



(26) 



In this approximation A''^ = A'^ = A and thus Gi — 0, 
so for the thermal transport lifetime we find, 

' . + ^2A 3m[,^H--pa/|0,^l)] ^^ 

2th \vf\ Jm 

(27) 

which agrees with results in Ref. |32| . We, however, aim 
to include the higher components of the order parameter 
expansion for a fully self-consistent calculation and for 
comparison with experiment. 



B. General properties of the thermal conductivity 
tensor 

As in I we focus on a tetragonal system with an 
open, along the z-axis, Fermi surface, and the mag- 
netic field applied in the xy plane, at angle 0o to the 
X-axis. We consider both d,j.2_y2 and dxy order parame- 
ters, and model their variation around the Fermi surface 
by y{(j)) = \/2cos2(/) and y{4>) = V^sm2(/) respectively, 
where is the angle between the projection of the Fermi 
momentum on the basal plane and the x axis. As before, 
we will consider both a cylindrical (no energy dispersion 
along z), and quasi-cylindrical (tight-binding dispersion 
along z) Fermi surfaces. The following considerations are 
valid irrespective of the Fermi surface shape. 

Experimentally, the in-plane (interplane) heat conduc- 
tivity is measured by driving the heat current along a 
high symmetry crystalline direction, such as [100] or [110] 
([001]). The longitudinal and/or transverse thermal gra- 
dient are defined and measured with respect to the di- 
rection of the heat current. This creates two physically 
distinct cases for the in-plane transport: the heat flow 
in the experiment is along either a node or antinode, see 
Fig.[2l If our task is, for example, to determine the shape 
of the gap from the measured thermal conductivity along 
the X-axis we cannot a priori assume whether the heat 
current is along a node or a gap maximum. 




FIG. 2: (Color online) Two distinct experimentally relevant 
setups for the thermal conductivity. Left panel: d^2_y2 gap 
symmetry, right panel: d^y gap symmetry. Response of the 
dxy superconductor to the thermal gradient along the [100] (x) 
direction with the field at an angle (pQ to this axis is equivalent 
to the response of a dx2_y2 system to the thermal gradient 
along the [110] (x') direction and the field at the angle (f>o = 
(j)' + n/4: to the x-axis. Note that the expriment is done with 
the heat current jh along [100], while the calculations are for 
the thermal gradient along this direction, see text for details. 



From the theoretical perspective, the knowledge of the 
full thermal conductivity tensor, Eq. (|23ap allows to de- 
termine the heat transport along an arbitrary direction. 
The two cases, dx2-y2 and dxy, see Fig.[2l transform into 
one another by rotation of the heat current: thermal con- 
ductivity measured for the heat current in the [100] di- 
rection for the dxy gap is equal to the thermal conduc- 
tivity for the dx2-y2 order parameter with the heat cur- 
rent in [110] direction. Therefore we focus on the d^2_j,2, 
and compute all the components Kij in the plane; the 
thermal conductivity for the dxy case is computed using 
these results. In a tetragonal system in the absence of 
the field, the off-diagonal elements Kxy = Kyx = 0, and 
the diagonal elements are equal, Kxx = Kyy, so that the 
conductivity is isotropic. Applying a magnetic field in 
the plane changes the situation dramatically. First, for 
the field applied at an angle 4>o relative to the [100] direc- 
tion, Kxx 7^ i^yy', it is casy to see (and we show it formally 
below) that Kyy{4>o) = Kxx(t^/2 + 4>o) since for these com- 
ponents the angle between the field and the heat current 
is the same. 

Second, the in-plane field breaks the tetragonal sym- 
metry, and therefore Kxy{4>o) 7^ for a general orientation 
of the field. We emphasize that this occurs even when the 
Lorenz force is neglected. The non- vanishing 

Kxy arises 

because of the difference between the transport along the 
vortex and normal to it; when the field (and the vortices) 
are at an arbitrary angle to the direction of the heat cur- 
rent, a transverse temperature gradient appears similar 
to the Hall conductivity in a material with the electric 
field applied at an arbitrary angle to inequivalent princi- 
pal axes. The transverse heat conductivity is of the same 
order or magnitude as the anisotropy between transport 
parallel and normal to the vortices, and hence much 
greater than the typical contribution proportional to the 
cyclotron frequency, oJc = eH/mc ^ {A'^ / E f){H / Hc2)- 
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Therefore the Hall angle is moderately large, see below. 

To make the argument more rigorous consider the gen- 
eral form of the heat conductivity tensor as a function of 
the field orientation, H. We write Eq. (|23ap in the form, 

i^viM - {vf4P>fAp) mpf ■ H)^ |A(p)p;T,i?)) . 

(28) 

Here the kernel /C((p/ • H)^, |A(p)|; T,i7) is determined 
by the equilibrium Green's functions, and, at a given 
point p/ at the Fermi surface, depends on the angle be- 
tween the Fermi velocity and the field, (p/ • H), the gap 
amplitude for that direction, |A(p)|, as well as on T and 
H. Since the kernel does not change if the direction of 
the field is reversed, we explicitly write it as dependent 
on(p/.H)2 

Let us start by considering a d^^-y^ gap. The inversion 
of the field in the xj/-plane corresponds to the change 
00 — *■ 00 + We can simultaneously change the variables 
in Eq. ((28|) according to {px,Py) — > {—px,—Py), which 
leaves the kernel invariant, and find 



Klji4>0 



nj((P0 



(29) 



for i, j = {x, y}; at the same time Kiz{(f>Q+TT) = —Kizifj^o)! 
and Kzz(0o + tt) = Kzz{4'o). Similarly, reflection of the 
field in the xz-plane, (j)o ~^ ~4'0i together with reflection 
Py —py again does not change the kernel, and leads to 



Hxxi-fjyo) 

Hyy{-<Po) 
Kxy{~4>o) 



Kxx{4>o) , 
Kyy(0o) , 
-Kxy{(t)o) . 



(30) 



Finally, if we rotate the field and the coordinate system, 
00 ^ 00 + Ti-ZZ [iPx,Py) {Py, -Px)]; wc find 



Kxx{(l>o + 7r/2) 
Kyy(0o + 7r/2) 

Ka;y(0O + 7r/2) 



Kyy(0o) , 
l^xx{4>0) I 
— Kxyi4>o) ■ 



(31) 



We carry out Fourier expansion based on these sym- 
metries. From Eq. Eq. and Eq. ^ we find 
for the d^2_,,2 gap 



^x~~y~ 

Kxxi4>o) 
Kj/y(0o) 
Kxy(0o) 



'So + ^^2 COS 200 + f^i COS 400 

Kq — K2 COS 200 -I- K4 COS 400 

^2 sin 200 -l- Kg sin 60o -|- . . . 



, (32) 



The cos 200 term in the longitudinal conductivity de- 
scribes the anisotropy between the transport along and 
normal to the vortices. Furthermore, if the superconduct- 
ing gap is isotropic (or absent), and hence the only depen- 
dence of the kernel K on the field orientation is via the 
term (p/ • H)^, it immediately follows that for our cylin- 
drically symmetric Fermi surface Hxx+^yy is independent 
of the field orientation, 0o, which requires K4 = 0. For an 
anisotropic Fermi surface there may be an angular mod- 
ulation of the thermal conductivity, but it would occur 
already in the normal state. Therefore the cos40o com- 
ponent that appears only in the superconducting state is 



predominantly due to the gap anisotropy. Such a decom- 
position in the analysis of the experimentally measured 
thermal conductivity was used in Refs. [l0lfl2ll4l to infer 
the gap structure of heavy fermion and organic quasi- 
two-dimensional materials; and in the following section 
we compare our results with their analysis. 

The origin of the sin 20o directional dependence of the 
transverse thermal conductivity is also transparent. In 
the presence of the field, the principal axes of the thermal 
conductivity tensor are along and normal to H. Conse- 
quently, when the heat current is along one of those axes, 
no transverse signal is generated, irrespective of the nodal 
structure. This is precisely the result found in high-Tc 
superconductors by Ocaiia and Esquinazi^^*^, who ob- 
served a nearly perfect sinusoidal thermal Hall response. 

We are now in the position to consider the differences 
between the d^2_y2 and dxy gaps. The longitudinal and 
transverse conductivities for the d^y order parameter are 
identical to the components of the thermal conductiv- 
ity tensor for the dx2-y2 case in the coordinates {x',y') 
rotated by a 7r/4 with respect to (x, y), see Fig. [2l 



'^x'y' 



(33) 



Moreover, the field applied at an angle 0o to the x' axis 
makes angle 0o = 0' + 7r/4 with the x axis, so that 
'^'i^'o) = 7?.(a)K(0o + 7r/4)7^~-^(a), with the rotation ma- 
trix 



n 



cos a sm a 
— sin a cos a 



(34) 



This leads to 

l^x'x' (0o) 
Ky'y' (0o) 
^^x'y' (0o) 



Kq + K2 COS 20Q 
Kq — K2 COS 20Q 

-K2 sin 200 -I- . . 



K4 COS 400 
K4 COS 400 



,(35) 



Importantly, this result implies that the four-fold term 
in the longitudinal thermal conductivity depends only 
on the orientation of the field with respect to the nodes 
of the gap. Indeed, let us restore the dependence on the 
angle 0o measured to the gap maximum, then 

Kx'x'iM = Ko + K2sin20o + K4COS40O + . . . , (36) 

Kx'y'iM = K2 COS 200 -I- . . . . (37) 

The last term in the longitudinal thermal conductiv- 
ity for the dxy order is identical to that for the dx2-y2 
gap. In other words, independently of the gap symme- 
try, the fourfold term in the longitudinal thermal conduc- 
tivity simply depends on the angle between the direction 
of the in-plane field, and the antinodal direction of the 
gap. Consequently, in the following sections we will fo- 
cus both on the overall features of the thermal transport 
and specifically on that term. 
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C. Calculated vs. measured thermal conductivity for the thermal conductivity, Eans. (|25p . (|26p . and ((271 



We will see below that the field-induced anisotropy in 
the transport along and normal to the vortices leads to 
the large thermal Hall angle, Kxy/i^xx ~ 0.1. In this 
case it is important to keep in mind that theoretical 
calculations are done under assumptions different from 
the typical steady-state experimental setup. The ther- 
mal conductivity tensor is defined via j/i..; = — k^ Vj-T, 
where jh is the heat current. The experiments arc done 
by driving the thermal current along a given ([100]) axis, 
while thermally insulating the sample in the transverse 
direction. The experiment measures the thermal gra- 
dients established under the conditions ]h,x = jh and 
hi,y — 0. Consequently, the measured longitudinal, k/ 
and the transverse, kj, thermal conductivities are 



VyT 

VxT 



Hi 



(38) 
(39) 



Presence of the off-diagonal terms does not substantially 
modify the absolute value of the longitudinal or trans- 
verse conductivity since K^y/ {KyyK^x) ~ 0.01 at most, 
and therefore ki ~ Kxx and Kt ~ Kxy 

Note, however, that our principal interest is in the 
fourfold nodal term, K4, which is itself only a frac- 
tion of the longitudinal thermal conductivity. Assuming 
Ho ^ 1^2, we find 



,(4) 



K4 



i%2 

2ko 



(40) 



In some region of the phase diagram, where K2 ^ k^, the 
two terms may be comparable. Our results indicate this 
range to be rather small. We find that the magnitude of 
the fourfold term is slightly changed by accounting for 
the difference between the computed and the measured 
quantity; however, the main features remain unmodified. 
Hence in the following we discuss the overall features of 
the thermal conductivity profiles, and only briefiy return 
to the difference between the computed and measured 
anisotropy in the conclusions. 



IV. CYLINDRICAL FERMI SURFACE 

Once again we begin by considering the anisotropy 
of the longitudinal heat conductivity, Kxxi4'o)i for a 
cylindrical Fermi surface, Vf — {vf cos <j>,Vf sin (j>,0) for 
< ^ < 27r and — tt/c < < tt/c, where c is the c- 
axis lattice spacing. As described in Ref. I, this FS does 
not allow for the self-consistent calculation of the order 
parameter in the vortex state for the in-plane field. As 
before, we restrict ourselves to the lowest order Landau 
wave function for the order parameter, take A{T,H) ~ 
A{T,0)[1 — H/Hc2\^^'^ , and use the corresponding results 



We choose to be direction independent, and carry 
out the self-consistent calculation in temperature and im- 
purity scattering according to Eqs. ipS]) and (P7|) . The im- 
purity self energy is determined in the unitarity limit with 
the normal state mean free path in/^o ~ 70, where ^0 is 
the coherence length. This toy model lends itself easily 
to numerical and, in some limits, semi-analytical work, 
and therefore allows investigation of the salient features 
of the behavior of Kxx- We show below that this model 
gives qualitatively correct results for quasi-two dimen- 
sional systems. In the self-consistent calculation there is 
a node-antinode anisotropy in the upper critical field at 
low T, and the comparison (given below) between the 
cylindrical and corrugated Fermi surfaces elucidates the 
role of this anisotropy for the behavior of the thermal 
conductivity. 

In general, to determine Kxx we need to self- 
consistently determine the DOS and the single particle 
lifetime as in the calculation of the specific heat, and 
then determine the transport time and the thermal con- 
ductivity. For the lowest Landau approximation these 
are given by Eqs. (|25II27[) . At finite energies this proce- 
dure can only be carried out numerically, as shown below. 
First we make some analytical estimates at low temper- 
atures, T = 0, and therefore set e = 0. We consider a 
dx^-y^ gap and focus on three values of the thermal con- 
ductivity: along the field, Kxxiff'o =0°), normal to the 
field, Kxxift'o = 90°), and for the field along the node, 
Kxx{(l)o = 45°). 

Define the mean free path £ = w/tq, where {2tq)^^ = 
—3mT,^{e = 0) is the single particle lifetime, which de- 
pends on the net density of states, and hence is sensitive 
to the direction of the field, 4>o, via the self-consistent 
T-matrix equation. The argument of the W-function 
and its derivatives in Eqs. (|25II27|) is then z — iA/£±, 
where £± = £\ sin((/) — <^o)|j and depends on the position, 
(j), on the Fermi surface. Since we work in the regime 
£ = VfTQ ^ A, we can set z ~ for most values of 
(j>, except for the directions nearly parallel to the field, 
\4> - 4>o\ < A/£ and |0 - 0o - 7r| < A/£. Let us denote 
the contribution from this narrow range as ki, and the 
contribution from the angles outside of this range as K2, 

so that Kxx = Ki + K2- 

We now estimate each contribution. For |(/) — 0o| < 
A/£ we use the expansion for large argument, W{z) sa 
i/{y/Trz) and W'{z) « —{/{y/irz'^) to estimate the contri- 
bution to the thermal conductivity as 



^i(0o) 
TNfvj 



4^' N Z'^"^''^' cos2 0d0 , 
3 J<j,o-A/e [1 + (Ato)2]3/2 

~"°^^"^ [l + (Aro)2]3/2 ' (42) 



where ai(0°) = 2A/^, ai(45°) = A/£, and ai(90°) = 
2A3/3£3. 

Over the remainder of the Fermi surface we set z ~ 
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and use W{0) = 1, W'{0) = 2i/^ to find for T ^ 



K'2{M 27r^ 



(43) 



cos 



27r 



= 8(AA/w/)2, (44) 
^^V^^^.^, (45) 



Here the prime denotes tliat we are integrating over the 
entire Fermi surface excluding the regions close to the 
field direction considered above. Notice that 3> 5^, 
and therefore the behavior of the thermal conductivity is 
controlled to much greater extent by the transport life- 
time than by the density of states. The transport lifetime 
is peaked along the nodal directions. 

These observations enable some analytical progress 
starting from the high fields, H < Hc2- In that case 
(5^ ^ /X <g; 1, and we approximate the density of states 
by its normal state value, Nf. Consequently, tq = t„, 
where r„ is the normal state scattering rate which has no 
dependence on the direction of the magnetic field. Defin- 
ing Kn = Tr'^Nfv'jTnT/3 we find angular variation of the 
thermal conductivity is approximately given by 



2cos2 



2tt 1 



cos^ 20 
I sin(0-0o)| 



(46) 



We assume here that fj, A/£, which is satisfied nearly 
everywhere up to Hc2 for clean systems. Now consider 
the behavior of the thermal conductivity for different di- 
rections of the field just below the upper critical field. 
For small /i we find the conductivity along the field, 
Kx2;(0°) ~ K„[l + (4/7r)/iln(u], while the conductivity nor- 
mal to the field is Kj,j;(90°) « k„(1 — 28/x/157r). Finally, 
for the field along a node, K^xi'^^") ~ k„(1 — 16/i/37r). 
Therefore in the immediate vicinity of the transition 
at low temperature we expect Kxx{QO°) > Kxx{4:5°) > 
KxxiO°), or nearly twofold profile of the thermal conduc- 
tivity. 

At lower fields and T = we enter the regime < 
1 <^ fi, where we can still approximate the density 
of states by the normal state value, but the thermal 
transport is restricted by sharp peaks in the lifetime for 
nodal quasiparticles, see Eq. (|46p . Linearizing the gap 
around the nodal points and carrying out the integra- 
tion, we find KxxiO°) ~ Kxxi90°) ~ Kn/i2^^^lJ-^/^), and 



Kxx{i5°) ~ i^n/i^n^^^)- Consequently in this regime we 
find Kxx{QO°) > Kxx{0°) > Kxx(45°), suggesting a weak 
minimum for the field along the node. Remarkably, at 
low T the conductivity normal to the vortex is always 
higher than that parallel to the field, but the amplitude 
of this anisotropy, and the relative position of the value 
of the thermal conductivity for the field applied along a 
node both change between /i ^ 1 and /i ^ 1. 
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FIG. 3: (Color online) Left panel: the Tif-diagrani of the 
longitudinal heat conductivity anisotropy, Hixx{4'o)- The ther- 
mal gradient is along the s-axis (maximal gap). Regions of 
different anisotropy are marked as [4>o,4'0)4'o]j which denotes 
Hxxi(t>o) > Kxxi(t>o) > Kxx{<f>o)- For points marked by circles 
and squares the profiles of the angle-dependent k,xx{4>o) are 
shown in Fig. [J] 



This analysis is supported by the numerical results. 
We show results only for the longitudinal thermal con- 
ductivity, da;2_j,2 gap, and the heat current along the 
antinodal direction. The complete heat conductivity ten- 
sor is given and discussed below for the corrugated Fermi 
surface, with self-consistent calculations; The results for 
both FS are very similar. 

The phase diagram of Fig. [3] shows regions with dif- 
ferent anisotropy of the heat conductivity. The changes 
along the vertical axis, T = 0, as a function of the field are 
in agreement with our estimates above: a.i H > 0.85i/c2 
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FIG. 4: (Color online) Center panel: evolution of the heat 
conductivity with temperature for H/Hc2 ~ 0.14. Right 
panel: heat conductivity as a function of the field angle for 
difi'erent fields at T/Tc = 0.25. The curves are shifted ver- 
tically for clarity. The framed numbers give true values of 
normalized Kxx{<i>o ~ 0) for several of them. Corresponding 
temperatures and fields are shown in each panel in the same 
order as the curves. 
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we find Kx2;(90°) > Kxx{4:5°) > Kxx{0°); while below 
that field Kxx{90°) > Kxx{0°) > K^a;(45°). Note that 
for H/Hc2 = 0.85 wc have 5^ ^ 1 and therefore already 
/i 3> 1. The variation of the thermal conductivity with 
the the direction of the applied field with respect to the x- 
axis (angle 0o) are shown in Fig.|4]for the points in the T- 
H plane marked by circles and squares in Fig. [3) Evolu- 
tion of with temperature at low fields (circles) is con- 
siderable: minimum for the nodal direction, (pQ = 45°, at 
low T quickly evolves into a maximum, and at T/T^ > 0.4 
the conductivity is largely twofold with no clear signa- 
ture of the nodal structure of the gap. The change of the 
anisotropy with the field at moderate T is more gradual 
(squares), and a pronounced peak at (f>o = 45°, for the 
field along the nodes, persists to moderately high fields, 
providing a clear signature of the nodal structure. 

The anisotropy between thermal transport normal to 
and parallel to the vortices is reversed at moderate tem- 
peratures (solid blue line in Fig. [3]): at low temperature 
i^xx{0) < Ka;a;(90) (in the notation of Eq. p2l) this means 
K2 < 0), while at high T wc find Kxx{0) > Kxx{QO) (or 
K2 > 0). This evolution is in agreement with that for 
a conventional superconductor, found by Maki2^. Note 
also that the four-fold term ^ K4 in Eq. p2l) is most 
pronounced at intermediate to low T and H. 



V. QUASI-CYLINDRICAL FERMI SURFACE 

A. Main results 

To solve the quasiclassical equations self-consistently 
for the order parameter, we need a model which allows 
for the c-axis superconducting currents when the field is 
applied in the a-b plane. Hence we analyze a corrugated 
quasi-2D Fermi surface given by 

P^^pI+pI- [r'^p)) cos{2s p,/r^pf) , 

so that the quasiparticle velocity has a nonvanishing 
^-component. This Fermi surface, with s = r = 0.5, 
was considered in I for the analysis of the specific heat, 
and we take the same values of parameters to directly 
compare the anisotropy of the heat capacity with that 
of thermal conductivity. Note that for this choice the 
DOS anisotropy in the normal state is ?^/(p) = 1, and 
the normal state conductivity anisotropy is / = 
s'^ = 0.25. For this anisotropy the vortex lattice is still 
Abrikosov-like. 

For the self-consistent calculation of the order parame- 
ter and Hc2 wc limit ourselves to three Landau level com- 
ponents in Eg. ljlOap . Ao,A2,A4, which is sufficient for 
the convergence of the calculation to high precision. As in 
I, we take impurity scattering in the unitarity limit with 
the strength in the normal state r/27r Tc = 0.007 (sup- 
pression of the transition temperature to Tc/Tco ~ 0.95, 
and the mean free path £n/^o — 70). We showed in I that 
this choice gives the following values of the critical fields 




FIG. 5: (Color online) Longitudinal heat conductivity for 
H||Vr||3;-antinodal direction: as a function of temperature 
for several fields (left); right: as a function of the field for 
temperatures {0.05, 0.2, 0.4, 0.6}. 



at T 



Tjantinode 



1.27So and 



~ 0.57^0, where Bq = ^'o/Stt^q where $0 = hc/2\e\. 
For the in-plane anisotropy we have then (^jj^nunode _ 
jjnodey jjanunode ^ ^^ie ratio between the c- 



axis and antinodal directions is 7J^2/-f^c2"**"°'^'^ = 0.4, 
similar to that observed in CeCoIns experimentally. 

We start by showing the temperature and field depen- 
dence of the heat conductivity tensor for the d^2_y2 gap. 
The heat current and the field are along the x-axis, along 
the gap maximum, as schematically shown in Fig. [3l 
The longitudinal thermal conductivity is seen in Fig. [S] 
to rapidly decrease below Tc{H) (left panel), as the gap 
opens in the single-particle spectrum. Notice, however, 
that the lines plotted for different fields intersect, imply- 
ing that the field dependence of Kxx is non-monotonic, 
as shown in the right panel. Kxx{H) increases with field 
at the lowest T and H. In this regime the low energy 
quasiparticles are located near the gap nodes, where or- 
der parameter vanishes, 3^(0„) = 0, and the transport 
lifetime, Eg. ipT)) is limited only by the impurity scat- 
tering, — UmE^. Hence in the competition between the 
increased number of heat-carrying quasiparticles due to 
field and scattering on the vortices, the density of states 
wins, and the conductivity increases with field. In con- 
trast, at higher T, the unpaired quasiparticles are already 
induced by temperature away from the nodes, and turn- 
ing on the field leads to increased scattering, hence the 
decrease in the thermal conductivity. The evolution of 
Kxx with T and H is nearly identical to that found for 
the field normal to the layers in a vortex state model 
with a single Landau level^, and is in agreement with 
experimental results^^. 

Fig. [6] shows the temperature (left panel) and field 
(right panel) dependence of the transverse heat conduc- 
tivity. As we emphasized above, for a general orientation 
of the field the dominant contribution to the transverse 
heat conductivity is not caused by bending of the quasi- 
particle trajectories in the applied field due to Lorenz 
force, but is due instead to the anisotropic scattering 
of quasiparticles by the vortices. If the thermal gradi- 
ent is not along one of these two "transport axes" , along 
the vortices and perpendicular to the vortices, transverse 
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FIG. 6: (Color online) Transverse heat conductivity for 
Vrj li-antinodal direction: H||noc?e as a function of temper- 
ature (left) and field (right). 



current arises. On the other hand, if the thermal gradi- 
ent is along or normal to the field direction, K^y = by 
symmetry, as is clear from Eq. ([5D|) . Hence vife show the 
transverse conductivity for 0o 45° which corresponds 
to H||no(ie. 

The transverse conductivity is allowed to change sign, 
as it only reflects the difference between the thermal con- 
ductivities parallel and perpendicular to H, which them- 
selves depend on the temperature and field. The tem- 
perature dependence shows a large peak at intermediate 
T for all H and tends to zero as the normal state is 
approached. The field dependence is more interesting. 
At low temperature, n^y is negative and monotonically 
decreasing up to the critical field, and rapidly goes to 
zero at lic2- At higher temperature K^y is positive, and 
has a peak at low fields. The temperature, at which the 
peak first appears, seems correlated with that where the 
downturn in K^xiH) is first seen, and the field value at 
the peak position moves in step with the minimum of 
Hxx{H) in Fig. [5] (right). It is therefore likely that this 
feature is a signature of the increased scattering due to 
magnetic field. This is supported by correlations between 
the peak and significant k(90°) — k(0°), which stems from 
magnetic scattering. 

To make connection to experiment, in Fig. [7] we show 
the temperature scans of the longitudinal heat conductiv- 
ity as a function of the field direction for low and mod- 
erate fields (left and middle panel respectively), and a 
field scan at T/Tc = 0.25 (right panel). The evolution 
with temperature at low field is similar to that found for 
the cylindrical FS, see Fig. [S] The low temperature re- 
gion is dominated by the evolution of the four-fold term, 
while, as temperature increases, the two-fold component 
becomes more prominent. The field scan strongly resem- 
bles the analogous result for the cylindrical FS: note the 
appearance of a pronounced peak for the nodal direction 
(00 = 45°) with increasing field. This shape of Hxx{<t>o) 
strongly resembles the experimentally found anisotropy 
in CcCoIus as shown in Ref[l(l. This speaks in favor of 
the d^.2_y2 gap symmetry in this material, and we will 
provide the detailed analysis at the end of this section. 

For the same relative orientation of the gap and the 




FIG. 7: (Color online) Heat conductivity profiles for T-scans 
at two different fields and a H-scan at T/Tc — 0.25. Dif- 
ferent symbols correspond to points in Fig. [10] (left). The 
curves are shifted vertically to fit on the same plot, and the 
value for 00 = for each of the shifted curves is given in the 
corresponding boxes. 



0.02 
0.015 
S 0.01 
0.005 


-0.005 



_T/T^ = 0.25 






0.065 




0.20 ■ 




0.26 




0.32 




0.39 




0.52 




0.65 













45 






0.05 




0.15 




0.20 




0.25 




0.30 




0.40 




0.6 




0.7 






FIG. 8: (Color online) Anisotropy profile of the transverse 
heat conductivity. This component gives the heat current 
along the {/-direction, when the VT||2: and the field makes 
angle 00 with the antinodal x-axis. 



heat current, we show a typical profile of the transverse 
thermal conductivity, Kxy, in Fig- [SI For 0o ~ 0° and 
00 = 90° this component vanishes identically, see discus- 
sion above. Over a wide range of T and H parameter 
range this component shows essentially sin 20o-like be- 
havior Eq. (|32p that agrees with experimental findings in 
high-Tc materials We emphasize that this modula- 
tion is completely unrelated to the nodal structure of 
the gap, moreover, only the deviation from the pure si- 
nusoidal profile, seen in several curves in Fig. [51 carries 
information about the gap structure. 

To complete the description of the thermal conductiv- 
ity, we present in Fig. [9] the longitudinal thermal conduc- 
tivity Kxx{4>o) for dxy symmetry of the order parameter 
(when the temperature gradient is along a nodal direc- 
tion), and compare it with the results for the rfj.2_y2 gap 
in Fig.[7l The temperature scans for field H/ Hc2 = 0.065 
in Fig. [7^a) and Fig. [5]^left) demonstrate that at high 
temperatures the same two-fold symmetry holds for both 
gap symmetries. At low T the fourfold feature disap- 
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FIG. 9: (Color online) Longitudinal heat conductivity 
anisotropy for d^y symmetry. The curves are shifted verti- 
cally, so the scale on the vertical axis shows only the absolute 
anisotropy. The relative anisotropy can be deduced from the 
true values of Kxx{0) shown in boxes. On the left is shown 
a temperature scan for H/Hc2 = 0.065, and on the right is 
a field scan for T/Tc — 0.25, shown correspondingly by filled 
circles and squares in Fig. [10] (right) . 



pears slightly faster for the d^y gap, and, at the lowest 
T, clearly has the opposite sign for the two gaps; this is in 
agreement with Eq. ([5^ and Eq. ([571) • Comparison of the 
field scans for T/T^ = 0.25, Fig. [Tjc) and Fig. [fright), 
shows a rather dramatic difference between the profile of 
the thermal conductivity under the rotated field for the 
two cases. A local maximum for the field at <j>o = 45° to 
the heat current is clearly resolved for dx2_y2 gap. For 
the same field direction either a minimum or no feature 
is seen for the d^y symmetry. Recall that at this tem- 
perature, T <C Tc, the thermally excited quasiparticles 
are still located near the nodes of the gap, at 0„. Their 
scattering on the vortices depends on the component of 
the Fermi velocity normal to the field, and therefore on 
the sine of the angle between the nodal and the field di- 
rections. Hence the strongest variation in the scattering 
occurs as the field sweeps through the nodal direction, 
when |0n — <^o| 1- For the dx2-y2 gap, with nodes at 
<l)n ~ 45° -|-n90° to the direction of the heat current, this 
leads to a noticeable feature in the profile of the thermal 
conductivity for 0o — For the dxy order parameter, 
with (j)n = n90°, the rotated field sweeps through the 
nodes at the same time as it is either parallel or nor- 
mal to the heat current, 4>q ~ 0°,90°. In that case the 
twofold transport anisotropy due to vortices masks the 
nodal signatures, and the signal is largely twofold. Only 
in restricted regions of the phase diagram, when the dom- 
inant twofold anisotropy is nearly absent (intermediate 
fields in the right panel of Fig. [H]), does the existence of 
the nodes affect the profile of k^x- This difference be- 
tween the behavior of Kxx{4'o) for the two types of gap 
strongly suggests to us that the experimental results for 
CeCoIns effectively rule out the dxy symmetry for this 
compound. 



FIG. 10: (Color online) The anisotropy phasediagram of the 
longitudinal heat conductivity k,xx{(I>o)- Notation [(j!>o, (jio, t/io] 
is short for hixx{4''o) > K-xx{(f>o) > Hxxi(t>o)- The solid blue 
line corresponds to the change of sign in the K2, which is 
positive at high T and negative at low T. Its position on the 
phase diagram is almost the same for the heat current along 
antinode and node. Left panel: the temperature gradient is 
along the antinodal direction. The symbols (circles, squares 
and triangles) correspond to anisotropy curves in Fig. [7] The 
variation of the transverse heat conductivity anisotropy for 
this gap orientation is shown in Fig. [8] Right panel: the 
temperature gradient is along a nodal direction. 



In Fig. [To] we summarize the results in the form of 
a phase diagram. The most noticeable differences with 
the cylindrical Fermi surface. Fig. [31 occur near the up- 
per critical field due to the Hc2 anisotropy, absent in the 
non-self-consistent calculation. Away from the critical 
field, in the low-to-moderate T,H corner, however, the 
anisotropy shows very similar features for the cylindrical 
and the corrugated FS, although the detailed positions of 
the separation lines, indicating the change in the shape 
of the thermal conductivity profile, is different. We be- 
lieve that the location of these lines is determined by the 
symmetry and the shape of the Fermi surface, and other 
microscopic details of the material. Recall that the cou- 
pling between different Landau level components of the 
order parameter is generated by the action of the differ- 
ential operator, v/(p) • V, which explicitly depends on 
the symmetries of the Fermi velocity. On the other hand, 
based on the similarities between the phase diagram com- 
puted with the lowest Landau level. Fig. [31 and that for 
three components, Fig. [TO] we conclude that the salient 
features and changes in the anisotropy as a function of 
temperature and field are captured here. 

It is also clear from the phase diagrams in Fig. [TUl 
and the anisotropy profiles in Figs. [7] and [9] that there 
is no simple relation between general shape and evo- 
lution of Kij(0o) for the two symmetries of the gap, 
dx2-y2 and dxy On the other hand, as we noted ear- 
lier, the coefficient K4 in the Fourier decomposition anal- 
ogous to Eq. ([3^ and suggested in Ref. [ij, KxxiM = 
f^o + K2 cos 200 -l- K4 cos 400 1 depends only on the orien- 
tation of the field with respect to the nodes. In Fig. [11] 
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FIG. 11: (Color online) The four- fold symmetry coefficient K4 
of heat conductivity for the scans shown above in figure [7] 




FIG. 12: (Color online) The T-H phase diagram for the four- 
fold component, K4, of the thermal conductivity under rotated 
magnetic field. Differently shaded regions correspond to dif- 
ferent sign of the four-fold term, and its connection to the 
nodal directions. The dashed line marks the crossover from 
< K,^(90°) at low T to K,,(0) > k^:.(90°) at high T. 



we plot the four-fold coefficient K4 for the vahies of T 
and H shown in Fig. [7] for the d.j;2_y2 gap. The four- fold 
anisotropy changes sign at intermediate temperatures for 
both small and moderate fields, see left two panels of 
Fig. [TlJ The coefficient K4 is small near Tc, where only 
the two-fold pattern that reflects the difference between 
transport parallel and normal to the applied field is de- 
tectable. The second sign change in the low-field range, 
left panel of Fig. [Til occurs at low temperature, close to 
the limit of validity of the BPT approximatior^i. How- 
ever, it is this feature that is connected in the phase dia- 
gram to the reversal shown in the middle panel for higher 
fields, see Fig. [121 which suggests that it is not an arti- 
fact of the approach, but a real effect. For the field scan 
at T/Tc = 0.25 (right panels of Figs. [7TTT|> K4 is always 
negative, so that the minima in the four-fold component 
mark the antinodal directions, while the maxima occur 
for H along the nodes. There is a sharp increase in the 
magnitude of this coefficient as we approach the critical 
field due to the in-plane Hc2 anisotropy. 



B. Comparison with experiment: CeCoIns and 
k-(BEDT-TTF)2Cu(NCS)2 

We close this section with a detailed comparison 
of our results with the experimental data. One of 
the main motivations for this work was the measure- 
ment of the thermal conductivity in a heavy fermion 
CeCoIns, — . Another example of quasi-2D supercon- 
ductor where the anisotropy was measured is k-(BEDT- 
TTF)2Cu(NCS)2-^. 

In CeCoIn5 the heat current is driven along [100] crys- 
tal direction. The observed profile of the heat conduc- 
tivity at T/Tc = 0.25 is in good agreement with that 
shown in Fig. [7] {dr^2_y2) for comparable temperature 
(T w 0.2Tc), including the peak for the field at 45° to 
the heat current. The profile differs significantly from 
that expected for a dxy gap as shown in Fig. [51 We 
find that the behavior of the experimentally determined 
four-fold term amplitude, K4, agrees with Fig. [TT] fright): 
it vanishes as ^ and saturates at H > 0.2 Hc2- 
The overall amplitude of this component is smaller in 
our computation than that observed experimentally by 
approximately a factor of three: however, since this mag- 
nitude is determined by the shape of the Fermi surface, 
we do not expect the model calculation to be quanti- 
tatively correct. Our results at moderate temperatures 
are consistent with the experimental temperature scan 
T = 0.15 - 0.9 Tc at « 0.1Hc2- At T - our results 
suggest an inversion of the K4-term which was not ob- 
served. However, in this region CeCoIns still has strong 
inelastic scattering (resulting in a peak of the thermal 
conductivity at T ^ 0.75Tc) which was not included in 
the calculation. Experimentally, extraction of the small 
K4 amplitude on the background of the dominant twofold 
term has greater relative errors. Moreover, since in this 
range kq 3> K2 ^ K4, the difference between the calcu- 
lated and the measured heat conductivity described in 
Sec. nil CI may also contribute to the discrepancy. Fi- 
nally, since the upper critical field in CeCoIns is param- 
agnetically limited, we cannot make a reliable connection 
of our results with experiment at low temperatures and 
high fields; in contrast, Zeeman splitting does not affect 
the low-to-intermediate T-H behavior. Therefore reli- 
able comparison can be made only in the region away 
from the critical field, where a maximum in the fourfold 
component points to the node, strongly implying c?^2_j,2 
symmetry in agreement with Ref . [lol. Note that, accord- 
ing to our analysis, generally the line of inversion of the 
fourfold term in Fig. [12] is distinct from the line separat- 
ing the increasing and decreasing k{H, T) at low fields, 
which was used to decide whether the minima or the 
maxima of the oscillations indicate the nodes in Ref. [lol . 

In the quasi-2D organic superconductor k-(BEDT- 
TTF)2Cu(NCS)2 the available data are in low-T, H re- 
gion onlyi^. The heat current is driven along [110] axis. 
Extensive analysis of the experimental data is required 
to separate the electronic contribution (which is small 
due low carrier density), from the phonon heat trans- 
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FIG. 13: (Color online) The T-H phase diagram for 
anisotropy of the heat capacity under rotated magnetic field. 
We choose for our model a quasi-cylindrical Fermi surface and 
a d-wave order parameter. For rotations of the field in the ah- 
plane the oscillations of the heat conductivity are four-fold, 
and a different feature (maximum or minimum) of this depen- 
dence determine the node position (arrows), depending on the 
location on this diagram. 



port. If, however, we concentrate our attention only on 
the behavior of the four-fold electronic term, the ob- 
served anisotropy fits well into the phase diagram for 
the dxy gap. Fig. [TOl In the region H ~ 0.07i?c2, 
T = 0.04- 0.07Tc, the fourfold term, K4COs4(/)o, has a 
maximum for 00 = 0° (along [110]) for T - 0.04Tc, and 
essentially disappears at T ^ 0.07Tc. In our mapping to 
the phase diagram, in the experimental regime the field 
along the node produces minima in the conductivity, and 
we concur with Ref. that the nodal direction is [100], 
which suggests the dxy symmetry. 



VI. CONCLUSIONS 

In this paper, following our approach in I for the spe- 
cific heat, and using the non-equilibrium Keldysh formu- 
lation of the quasiclassical theory, we derived a general 
expression for the heat conductivity tensor of a supercon- 
ductor in magnetic field. The derivation was based on the 
closed form solution for the Green's function obtained in 
I, that made use of the Brandt-Pesch-Tewordt approx- 
imation. The utility of this approach lies in its ability 
to self-consistently take into account impurity scatter- 
ing, the detailed shape of the Fermi surface, and multi- 
ple Landau levels in the order parameter in the vortex 
state. Numerical computations based on this approach 
are very time-efficient. The main advantage of our ap- 
proach is that it provides unified method for calculations 
of transport and thermodynamics over a large range of 
temperatures and fields, well beyond the realm of appli- 
cability of the semiclassical schemes. 

In these two papers we applied the developed formal- 
ism to a d-wave superconductor with a quasi-cylindrical 
Fermi surface. We concentrated on the behavior of the 



specific heat and thermal conductivity, since they are the 
most widely used experimental probes. To make a con- 
nection between theory and experiment, we provided, for 
the first time, a complete description of the anisotropy 
of the thermal conductivity and specific heat in the T-H 
phase diagram, starting from the "semiclassical" region 
at low r, H and up to the critical field. 

Two figures summarize the main results of this work. 
Figure [12] recalls the results of I, and shows the phase 
diagram for anisotropy of the specific heat under rotated 
magnetic field, C((/)o). The preceding figure. Fig. [T^ 
shows the anisotropy of the fourfold component, K4, of 
the longitudinal thermal conductivity, Kxx , for the model. 
One of our main findings is that both anisotropic signa- 
tures change sign, i.e. invert, in the T-H plane. For the 
specific heat, the inversion of the anisotropy is due to the 
effect of the quasiparticle scattering on vortices on the 
density of states, see I. The semiclassical (Dopplcr shift) 
picture predicting a minimum of C for the field along 
a node is valid at low T and H, where it was designed 
to work. The effects of the energy shift, however, are 
superseded by the redistribution of the spectral weight 
due to scattering, which becomes dominant not only at 
moderate fields, but also at low fields for finite energies. 
Consequently, the anisotropy changes sign at finite T, 
Fig.m 

Analysis of the heat conductivity is more involved due 
to interdependence of the transport scattering time and 
the density of states in the self-consistent treatment. We 
showed that, under a rotated field, the four-fold term in 
Fourier decomposition of the heat conductivity, «;4(0o)j 
exhibits signatures of the nodes, and depends only on 
the angle between the field, H, and the nodal directions, 
but not on the orientation of the heat current relative 
to the nodes. From comparison of Fig. [T^] and Fig. [T3] 
it is clear that the evolution of the fourfold coefficients 
in the specific heat and thermal conductivity, including 
sign changes, is quite similar across the phase diagram. 

The exact location of the inversion lines depends on 
the microscopic details, such as the Fermi surface shape 
and the detailed structure of the order parameter. The 
relative position of these lines, however, is stable with re- 
spect to the moderate changes of the FS curvature along 
z-axis and impurity concentration. The developed the- 
ory is valid over most of the phase diagram, except for 
very low fields in dirty samples, where the averaging pro- 
cedure in the Brandt-Pesch-Tewordt approximation is no 
longer valid. Thus the qualitative changes in the fourfold 
term at moderate fields, where the anisotropy is the most 
prominent, should be detectable experimentally (albeit 
may prove labor-intensive). One possible experimental 
approach to search for the node locations would be to 
measure anisotropy at several points in the phase dia- 
gram, in order to map out the evolution of the anisotropic 
contribution. 

Finally, by comparing our results with available experi- 
mental data on the specific heat and thermal conductivity 
we concluded that the order parameter of heavy fermion 
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material CcCoIng^iiS has d^'^-y'^ symmetry, and recon- 
ciled the thermodynamic and transport measurements. 
Analysis of the thermal conductivity for organic super- 
conductor K-(BEDT-TTF)2Cu(NCS)2i^ places it most 
likely into dxy family. We believe that our method, which 
allows detailed microscopic calculations for specific com- 
pounds, will enable unambiguous interpretation of the 
anisotropy in thermal and transport properties of uncon- 
ventional superconductors, and will lead to maturing of 
this method as a tool for determining the nodal direc- 
tions. 
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APPENDIX A: DETAILS OF THE HEAT 
CONDUCTIVITY DERIVATION 

Our approach allows to obtain an expression for ther- 
mal conductivity that generalizes previous results to the 
case of vortex state with multi-Landau level order pa- 
rameter, and arbitrary impurity strength. 

Keldysh diagram technique is formulated for 8x8 
Green's functiouf^^ which is traditionally split into three 
4x4 parts. Retarded (R), Advanced (A) and Keldysh 
(K):g^,X = {i?,A,X}, 







(Al) 



g (p,R,e)- (^^^^^^^^x^^ 5^+g^<T* 



(A2) 



For stationary problems these functions obey the normal 
ization conditions: 

r^^r-^^-^n , 5^3^ -f 5^ = 0. (A3) 

We do not need to solve equations for all the functions 
as they are related through symmetries,— 

5^(p; e) = g^iv; e)* g^(p; e) = g^(p; e)* 
3^(p; £) = 5^(-p; -e) g^(p; £) = g'^J-p; -e) 
5^(p; e) = 5^(-p; -e) g^(p; £) = g^'(-p; s) 

(A4) 

/^(p; e) = /^(p; ef f«(p; e) = f^(p; e)* 
/«(p;£) = /^(-p;-£) f«(p;£) = f^(-p;-£) (A5) 
/^(p;e) = -/^'(p;£)* f^'(P;^) = 



The different functions obey transport equations: 

[(£ - cjb)9z - + i^f^g^^^ = (A6) 

((e - (tb)% - - ?^((e - aB)% - a^) 

-a^g^ + g^?^ + ivfVg^ = . (A7) 

Here 

(A8) 



e . 
--Vf A , 

c 



is the coupling of quasiparticles to an external magnetic 
field. The self-energy is decomposed into the mean-filed 
order parameter and the impurity contributions, a'^ = 
Self-consistency equations for the singlet 



imp ■ 



order parameter are 

A^Xp,R) = A^(p,R) 



(A9) 



J ^ Jdp'nfip')V,ip,p')f^ip',R;e), 



A^(p,R)==0, 



(AlO) 



and for triplet superconductivity the equations are iden- 
tical, upon replacing by its vector counterpart, A"**^. 
The Keldysh part of the self-energy comes in this case 
from impurities only. The self-consistent t-matrix ap- 
proximation for isotropic impurity scattering gives 

dfrnpi'^; £) = n,mpt''(R; e) (All) 
P'^ = lii + uNf{g^-^)t^-'^ (A12) 
^ Nft'^ig^)?'^ (A13) 

where the angular brackets denote the normalized Fermi 
surface average as in the main text of the paper. 

In this appendix we denote functions in thermal equi- 
librium by the subscript 'eq', but in the main text we 
omit it for brevity. In local equilibrium, 

5^ = g^'^^. ^^.gf. , = ^?.1*o, - , (A14) 



Heat current is 



= tanh 



+ 00 



2T(R) 



(A15) 



hiR)^2Nf J dpfUfipf) J ^ev/(p/)5^(p/,R;£), 

— OO 

(A16) 

where g^(p/,R;e) is the diagonal part of 'g^ . Using 
the Green's function symmetries one can show that in 
this formula it is equivalent to g^ = l/4Tr(g^). The 
factor of two reflects our definition of Nf for single spin 
projection. Thermal conductivity in the linear response 
is determined from 

6h{R) = ~kVT^ (A17) 

4-00 

= 2Nf J dpfUfipf) J ^ev/(p/)lTr[,Sg^(p/,R;£)]. 
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Our goal is to find the linear, in the temperature 
gradient that drives the heat current, correction to the 
Keldysh Green's function, 5g^ . The Green's functions 
varies on the scale of the magnetic length (or intcrvortex 
distance). A, and on the scale of the superconducting co- 
herence length, ^0- We assume that A,<^o ^ Lt, where 
Lt is the length scale for temperature variation. We now 
write the gradient term as a sum of the gradients due to 
inhomogeneity and due to the externally imposed slow 
temperature variation. 



d_ 

df 



(A18) 



J 



with the last term much smaller than the first. 
Solution in the local equilibrium is obtained from 



R,A 



(A19) 



We note that the following analysis can be easily adjusted 
if we add an external potential to hf^-^. For now we 
continue without it and write down linearized equations 
for Green's functions near equilibrium, g-^ 
with driving term due to v/Vt = Ot, 



5^ +' 



[h. 



R,A 



K 



R,A\ 

g'^5d^ -Sd'^g^ 

■-'cq ■-' oq 



R,A 
K i 



JR^K 



We decouple equations for Sg^'^ from 57/^ by introduc- 
ing Eliashbcrg anomalous propagato r^^i^^ and the self- 
energy. 



^cq ^cq 



The heat current is determined by 6^° 



(A22) 
(A23) 



which satisfies 



(A24) 



= -zv; VT$oq (sfq - dfj + Sd^'gi - g^^d'^ . 

We need to solve this equation together with the self- 
consistency equations on . Normalization requires 



i/cq 



0, 
= 0. 



(A25) 
(A26) 



Up to this point all the equations are valid for both sin- 
glet and triplet pairing states. Below we focus on the 
singlet pairing and only briefly comment on the differ- 
ences between the singlet and the triplet cases. The self- 
consistency for the retarded and advanced linear correc- 



J 



ivfVxiSg^) + ivfVTg. 



K 



(A20) 
(A21) 



tions contain the order parameter and the impurity con- 
tributions, = 5A + Sa^^p, with 

5A = J^J dp'n/(p')F(p,p')<5/^(p',R;e)(A27) 

n.™p7V/ff/(<S5'''-'>?q^,(A28) 



R,A 
imp 



but the anomalous part is due to impurities only, 

5dl,p = n^mpNf f« {5r)ti (A29) 

We see that the equations for the anomalous Green's 
function, (5g°, and the self-energy, (5(t°, arc completely 
decoupled from those for the retarded and advanced 
Green's functions. On the other hand, the equations for 
ggR,A ^QpQYif^ on the anomalous 5g^ through the vari- 
ation SA. For simple retarded and advanced Green's 
hmction {Y = R,A), 



9c<n •/cq 

we obtain the impurity t-matrix in equilibrium. 



(A30) 



r sin^ So 



w 1 - ^((grq)^ - ifi) irj + 



r 



za, {(fl)/n) ctg^o - (.9.0/^ 



(A31) 



Here we solved (|A12p using the BPT approximation. vanish due to inversion symmetry, p —p. 

Note that for the triplet case the off-diagonal parts (/oq) , , , , , 

For the smglet case we assume oa — 0, and validate 
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this assumption at the end of the calculation. We show 
that the linear correction, 67)°', a product of an even, 
under p — p reflection, function and an odd, in mo- 
mentum, term v/(p)Vt, so that its average over the FS 
vanishes. 

We only need to solve the equation for the anoma- 
lous part of the Green's function, since neither Retarded 
nor Advanced part has a unit-diagonal term. Hence 
Tr[g^^"^] = 0, and they do not contribute to the heat 
current. Now Tr[(5g^] = Tr[(5p°], and we find 

= -*v/Vt$., (?f,-?^J (A32) 

where the matrix, 

oq V Gq J-* / iJ t!q cq 

= efl,^f3 _ nyi - , (A33) 



depends only on the equilibrium self-energies, (A^^ = 
for a triplet) 

iR^A^^^^R,A_^^^ (A34) 

^oq' = ^impt^' J (A35) 
^(!^'^ ~ T^imptS , (A36) 

Af;^ = A„q + n^rapt^^ . (A37) 
We parameterize the anomalous propagator as 

and find the equations for the diagonal components, g° 
and g"\ from (|X32)) 



-.9"i?" + - e^) + \r{^l - Afj + (A« - Af^) + ^v;V, g°^G, (A39) 

-g'°D° + - e^) + Afj + \r{K + K) + V,, 5'" = -»v; Vt <f>„, (sf, - 5.^J , (A40) 

where we defined D° = — D\ Expressions for the off-diagonal terms arc obtained from the normalization 
condition (|A26|) . 

r = -n^i-g'-iC + /4) + ff"(/e^ - /el)) (A41) 

5" cq 9 cq 

r = ^^^{9%f'' + /■") + - /^)) (A42) 

Q-Tt Q/i cq cq' cq cq ' ' 

cq i'cq 

Combining these equations, and using the BPT approximation (i.e. assuming that g^^g'^ =const and spatially 
averaging the terms containing /, / and A), we obtain the final expression for the unit-diagonal part of the anomalous 
propagator, 

9" = g^gZ^O^S-^ - 4)(*^/ V*eJ , (A43) 
with the following definitions of coefficients, 

1 





= -D" 


^ 2(.9. 


G2 


= e^- 


e^ + , 


G3 


= 




G4 


= -D" 





i-ifl!^ + /4)(A., - A J + (/ + OiK - K)] (A44) 



1 



[(/.^- - - A J + (/f^ - tj(^^^ - ^^J] (A45) 

^ [(/e^^ + /4)(Af, + Afj + (/f + /l)(Af, + A^J] (A46) 



2{9S^-9^) 



2{9S^-9i) 
1 



y [-(/.^^ - /4)(A., + A J + (/:: - r )(Af, + A^^J] . (A47) 

In order to prove that our assumption of Sa^ = is justified, we note that the cocfhcicnts Gi are even under inversion 
of p. Then the resulting g° and g'° arc odd due to additional factor vj V$^q, and their averages over the Fermi surface 
vanishes. The same is true for the off-diagonal functions /° and /° in singlet case. This assumption is not valid for 
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triplet pairing: even though g"" and g'°' are odd under inversion of p, wc sec from Eqs. (|A4ip and (jA42p that /° and 
/° are even (since /„q's are odd), and additional terms due to 5d°- appear, making the self-consistent solution of the 
equations more difficult. Exceptions to this statement exist for certain order parameters and for special orientation of 
the temperature gradient. For example, when VT is applied in a direction along which A does not change, we find 
(/°) = 0. Two obvious examples are a) ^{pz) and VT is in the xy-plane and b) A{px,Py) and VT || z. ) 
The resulting expression for the heat conductivity is 



T 



Nf 



de 



47r 7 T VT2cosh''(£/2T) 



Go 



(A48) 



r 



We checked that this expression for a uniform supercon- 
ductor agrees with the heat conductivity of Graf et al^. 
This completes the derivation of the heat current. We 



remind the readers that in the main text we drop the 
equilibrium subscript 'eq' to make expressions less clut- 
tered. 
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